/*  sdflogag.cpp  by K.Tsuru  */
// function ID 3405  DARDIX. Remade ver. 2.30
/****************************************************
log(x) by AGM (Arithmetic-Geometric Mean) method for x > 0
It needs the constant Pi().

[Reference]
"Arugorizum Jiten(Algorithm Dictionary, in Japanese)"
Kyoritsu Syuppan, 1994, p.363
*****************************************************/
#ifndef SN_H
#include "sn.h"
#endif
static void Kt3(const SDouble& q, SDouble& K, SDouble& t3){
    SDouble t0(q), qnn(q), q2(q*q), q2n1(q2*q);
    t3 = q;
    int i = 0;
    while(1){
        qnn *= q2n1;
        if(i & 1) t0 += qnn;
        else      t0 -= qnn;
        t3 += qnn;
        if(qnn.IsMLT(ONE)) break;
        q2n1 *= q2;
        i++;
    }
    t0 = ONE - 2*t0;    t3 = ONE + 2*t3;
    q2 = t0/t3;     q2 *= q2;   q2*=q2;
    t0 = ONE - q2;  K = Sqrt(t0);
}
/***********
x must be 0 < x < O(1)
************/
static SDouble SNLogAGM(const SDouble& x){
    RealSize C;
    SDouble a(1.0), b, q, d, t3, K, r;

    Kt3(x, b, t3);  //It evaluates 'b' and 't3' from 'x'.
    while(1){
        q = (a+b)/2.0;
        b = Sqrt(a*b);
        a = q;
        d = a-b;
        if(d.IsMLT(a)) break;
    }
    a = (a+b)/2.0;
    d = t3*t3;
    r = Pi()/(a*d);

    r.ChangeSign();
    return r;
}
/***********
main body
***********/
SDouble LogAGM(const SDouble& x){
    SDouble X, add;
// x = X*10^exp. log(x) = log(X)+exp*log(10), add = exp*log(10), 0< X < 1
    if( GetLogxCalcMethod(x, X, add) ) return X; // x= 1 or 10, X=0 or log(10)
    return ( SNLogAGM(X) + add );
}